%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  This is a paper in REVTeX, submitted to J. Appl. Phys., March.31, 1998.
%
%   14 figures:   1. flowchart.ps        2. noise.ps
%                 3. freq5.ps            4. corr5.ps
%                 5. per_corr_1700.ps    6. per_corr_760.ps
%                 7. per_corr_284.ps     8. signal.ps
%                 9. spectra1000.ps     10. corr1000.ps
%                11. spectra100.ps      12. corr100.ps
%                13. spectra100_c.ps    14. corr100_c.ps
%
%     in PostScript format at the end, following the .tex file.
%
%    The figures can be easily included in the printout by getting rid
%      of the comments before ``newcommands'' \psA, \psB and \PSbox
%
%  The entire manuscript in PostScript (with figures) can be accessed at 
%                 http://mmm.mit.edu/~liju99/cam.ps
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%     << Spectral Method in Thermal Conductivity Calculation >>
%
%                      Ju Li, Sidney Yip
%             Department of Nuclear Engineering
%    Massachusetts Institute of Technology, Cambridge, MA 02139
%
%                    email: syip@mit.edu
%                   phone: (617) 253-3809
%       mailing address: Room 24-208, MIT, Cambridge, MA 02139-4307
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% \documentstyle[eqsecnum,preprint,prb,aps]{revtex}
\documentstyle[eqsecnum,prb,aps]{revtex}
\def\btt#1{{\tt$\backslash$#1}}

\newcommand{\PSbox}[3]
{\begin{center}\mbox{\special{psfile=#1}\hspace{#2}\rule{0pt}{#3}}\end{center}}
\newcommand{\psA}[1]{\PSbox{#1 hoffset=-40 voffset=-15 hscale=45 vscale=45}
{7.5cm}{17cm}}
\newcommand{\psB}[1]{\PSbox{#1 hoffset=-98 voffset=-200 hscale=100 vscale=100}
{15cm}{18cm}}

\begin{document}
\draft
 
\title{Spectral Method in Thermal Conductivity Calculation}
\author{Ju Li, Sidney Yip}

\address{Department of Nuclear Engineering \\
Massachusetts Institute of Technology, Cambridge, MA 02139-4307}

\date{\today}
\maketitle
\begin{abstract}
  We discuss the numerical issues concerning molecular dynamics
  calculation of thermal conductivity using the Green-Kubo formula.
  For crystals with high thermal conductivity such as SiC, converting
  the MD heat current data into thermal conductivity is non-trivial.
  One can accelerate the process using Fast Fourier Transform and the
  spectral method, and we shall clarify some mathematical subtleties
  that may arise. Due to finite data length, integration of the
  correlation function must be terminated before noise takes over; we
  propose two termination criteria which give reasonable results for
  crystalline $\beta$-SiC.  In the end we derive a simple solvable
  model that illustrates the above points and can be used to check the
  codes.
\end{abstract}

\pacs{PACS number: 02.70.Ns, 33.15.Vb, 66.70.+f}
 
\narrowtext

\section{Introduction}
\label{sec:introduction}

We have calculated the thermal conductivity tensor
$\underline{\underline{\kappa}}$ of solid SiC with and without
defects\cite{Li} using the Green-Kubo formula\cite{Green-Kubo}
\begin{equation}
\underline{\underline{\kappa}} =
\frac{1}{k_{\rm B}T^2\Omega}\int_{0}^{\infty}d\tau \langle
{\bf J}^q(0){\bf J}^q(\tau) \rangle,
\label{first}
\end{equation}
where ${\bf J}^q$ is the instantaneous heat current of the system,
provided by a molecular dynamics simulation
run\cite{Ladd,Lee,Kitagawa,Kaburaki}. Details about the interatomic
potential model, simulation method and results for this material can
be found at the above reference. Our main concern in this paper is one
technical aspect which is often not fully treated by various
literatures, that is how one evaluates
$\underline{\underline{\kappa}}$ once the MD simulation finishes and a
finite length of the ${\bf J}^q(t)$ data has been recorded.

The problem is not critical for thermal conductivity calculations in
liquids and non-crystalline solids simply because the heat current
correlation function
\begin{equation}
\underline{\underline{g}} (\tau) = \frac{1}{k_{\rm B}T^2\Omega} \langle
{\bf J}^q(0){\bf J}^q(\tau) \rangle
\label{corfun}
\end{equation}
has a rather short correlation time, as the heat carriers quickly get
scattered. Thus it is a relatively quick calculation which does not
require very long simulation runs to get convergent results, and the
correlation function can even be evaluated on the fly.  However, a
perfect crystalline solid material usually has a much longer
correlation time (by one or two decades), and thus a much higher
thermal conductivity\cite{explain_g}.

In the case of crystalline $\beta$-SiC, even at the high temperature
of 1500K, to get a barely convergent result of the correlation
integral requires a MD run longer than 2 million timesteps (800 $ps$).
Further taking into the account that
\begin{enumerate}
\item To get realistic results, an interatomic potential of some
  sophistication, in our case the Tersoff three-body
  potential\cite{Tersoff}, must be used for this industrial material.
  
\item The calculation of thermal conductivity depends sensitively on
  the simulation cell size, due to the intrinsic mechanism of
  three-phonon scattering. In our work we find that a 216-particle
  cell is the minimum, while 512 or 1000 particle cells are
  recommended.
  
\item At low temperatures, longer simulation runs need to be taken
  because the correlation time gets even longer. For instance, at 500K
  a minimally convergent run is 8 million timesteps, or approximately
  4000 $ps$.
  
\item We find that the result of a simulation run is sensitively
  dependent on the {\em initial} particle velocity (phonon density)
  distribution of the run, because when the crystal has a high thermal
  conductivity the existent phonons are not easily scattered; the
  system will retain memory of the initial condition for a long time,
  staying in a limited region of phase space and the ergotic
  hypothesis not fully coming into effect at least for reasonable run
  lengths. For this reason, we need to do several independent runs
  starting from different initial conditions, and average over the
  respective correlation functions.
\end{enumerate}
the necessary computational effort is quite formidable. Kitagawa {\em
  et al}\cite{Kitagawa} have reported similar constraints in
calculating the thermal conductivity of AlN.  Because of this, the
treatment and analysis of the raw data, which sometimes is barely
sufficient, becomes an important issue.

The computation of $\underline{\underline{g}}(\tau)$ itself can be
quite expensive. For instance, in the above case at $T=1500$K and
$N=2\times 10^6$ timesteps ``minimal-run'', the correlation in
$\underline{\underline{g}}(\tau)$ persists for $\tau < 20$ $ps$, which
corresponds to $K=60,000$ timesteps. If one tries to directly compute
$\underline{\underline{g}}(\tau)$ by brute force averaging, he faces
multiplication operations on the order of $NK=10^{11}$, which means
evaluating all tensorial components will take more than 2 hours if the
subroutine runs at 100 Mflop/s.  One can of course down-sample the
data or decrease the number of time origins, but getting the raw data
is even a more expensive process that one certainly does not want to
waste any information here. Furthermore, one is not sure what the
correlation time is beforehand, and wants to see the entire
$\underline{\underline{g}}(\tau)$.

The spectral method delivers more and is faster, by taking advantage
of the Fast Fourier Transform algorithm\cite{Press}.  The number of
operations is ${\cal O}(N\log N)$ with a small prefactor on the order
of 10. Because $\log N$ is always much smaller than $K$, it is faster
than the brute force method. It does not waste any information
residing in the raw data; and it also provides frequency space
information.

A flowchart of the procedure is shown in Fig. \ref{flowchart}.
Basically one takes the FFT of the ${\bf J}^q(0<t<T_{f})$ data to get
${\bf J}^q(\omega)$, square it and divide by $T_{f}$ to get the power
spectrum $\underline{\underline{g}}(\omega)$, and do an inverse FFT to
get $\underline{\underline{g}}(\tau)$. One can either directly read
off the thermal conductivity by taking a limit of the power spectra
$$
\underline{\underline{\kappa}} = \lim_{\omega \rightarrow 0} 
\frac{1}{k_{\rm B}T^2\Omega} 
\frac{{\bf J}^q(\omega){\bf J}^{q*}(\omega)}{2T_{f}},
$$
or take the more cautious route of integrating
$\underline{\underline{g}}(\tau)$, whereas a termination of the
integration is needed sooner or later due to finite data length.  We
shall propose two termination criteria which works for solid SiC, that
minimizes human intervention, yet giving reasonable results.

The spectral analysis of continuous variables will be introduced in
Section \ref{sec:spectral}. Its discrete implementation will be
explained in Section \ref{sec:discrete}, and the termination criteria
for SiC will be explained in Section \ref{sec:termination}.  A simple,
analytically solvable model
$$a_i = a_{i-1}e^{-\frac{1}{L}}+ R_{i}$$
is proposed in Section
\ref{sec:model}, which illustrates most of the above issues and can be
used to test the code, hopefully providing a handle for beginners.

\section{Spectral Analysis}
\label{sec:spectral}

Consider a truncated function $f(t)$ which is finite from $t=0$ to
$T_{f}$, and $0$ elsewhere, so it is absolutely integrable. Then
\begin{eqnarray}
f(\omega) &&= \int^{+\infty}_{-\infty} f(t) \exp(i\omega t) dt,
\label{first:1} \\
f(t) && = \frac{1}{2\pi} \int^{+\infty}_{-\infty} f(\omega)
\exp(-i\omega t) d\omega.
\label{first:2} 
\end{eqnarray}
Notice that the artificial truncation of $f(t)$ will come to affect
$f(\omega)$ statistically only when $\omega T_{f}$ is comparable to
unity. In order to rule out the arbitrariness and make sure that the
information we get is intrinsic to the signal, we should take care not
to use the spectrum between $0<|\omega|<\pi /T_{f}$, which is of
course a very small region because $T_{f}$ is large, and will vanish
as $T\rightarrow \infty$. Nevertheless, we cannot use $f(\omega)$ at
$\omega$ equals absolute zero, because all our applications have
finite $T_{f}$. We can only use $\omega$'s which are much greater than
$\pi /T_{f}$ but are still very small, as the limiting behavior of
$w\rightarrow 0$.

The correlation function of $f(t)$ can be defined to be
\begin{equation}
g(\tau) = \frac{1}{T_{f}} \int^{+\infty}_{-\infty} f^*(t)f(t+\tau) dt.
\label{last}
\end{equation}

For $g(\tau)$ to behave well, there should be $\tau \ll T_{f}$. Let's
define a maximum correlation length $\tau_{\rm max}$ above which
$g(\tau)$ is negligible; thus $\tau_{\rm max} \ll T_{f}$.

From Eq. (\ref{first:2}), (\ref{last}), there is
\begin{eqnarray} 
g(\tau) =&& \frac{1}{4\pi^2 T_{f}} \int^{+\infty}_{-\infty} d\omega
\int^{+\infty}_{-\infty} d\omega^\prime f^*(\omega) f(\omega^\prime) 
\exp (-i\omega^\prime \tau) \nonumber\\
&& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
 \cdot \int^{+\infty}_{-\infty} \exp(i(\omega-\omega^\prime)t) dt 
\nonumber\\
=&& \frac{1}{4\pi^2 T_{f}} \int^{+\infty}_{-\infty} d\omega 
\int^{+\infty}_{-\infty} d\omega^\prime f^*(\omega) f(\omega^\prime) 
\exp (-i\omega^\prime \tau) \nonumber\\
&& \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;
 \cdot  2\pi{\rm \delta} (w-\omega^\prime) \nonumber\\
=&& \frac{1}{2\pi T_{f}} \int^{+\infty}_{-\infty} d\omega |f(\omega)|^2 
\exp (-i\omega \tau) 
\label{G}
\end{eqnarray}
which means that the Fourier spectrum of $g(\tau)$ is just the {\em
  power spectrum} of the signal $f(t)$.

There is a frequently quoted result that the integral of $g(\tau)$
from $0$ to $+\infty$ is just the power spectrum $|f(\omega)|^2/T_{f}$
taken at $\omega = 0$ and divided by $2$\cite{Lee}. This needs some
clarification, because if we put $\omega = 0$ in Eq.
(\ref{first:1}) directly, then
\begin{equation} 
f(\omega=0) = \int^{+\infty}_{-\infty} f(t) dt = \int_{0}^{T_{f}} f(t) dt,
\end{equation}
which is a random number of no physical significance.

This is not what we expect. To see what is wrong, let
\begin{equation} 
A = \int^{\tau_{\rm max}}_{-\tau_{\rm max}} g(\tau)  d\tau,
\end{equation}
then by Eq. (\ref{G}),
\begin{eqnarray} 
A &&= \frac{1}{2\pi T_{f}} \int^{+\infty}_{-\infty} d\omega |f(\omega)|^2
\int^{\tau_{\rm max}}_{-\tau_{\rm max}} \exp (-i\omega \tau) d\tau \nonumber\\ 
&&= \frac{1}{2\pi T_{f}} \int^{+\infty}_{-\infty}  
d\omega |f(\omega)|^2 
\frac{2\sin(\omega \tau_{\rm max})}{\omega}. 
\label{puff}
\end{eqnarray} 

Since $T_{f}$ can be very large as we do long simulation runs, $\tau$ can
also be allowed to be very large, {\em provided that} $\tau \ll T_{f}$.
Thus
\begin{equation} 
\frac{2\sin(\omega \tau_{\rm max})}{\omega} \rightarrow 
2\pi \rm{\delta}(\omega)
\label{now} \end{equation} 
as $\tau_{\rm max} \rightarrow +\infty$. Note that however, as we said
before, we cannot really use the {\em absolute} $\omega=0$ value of
$f(\omega)$, but only its limit well above $1/T_{f}$, so the system is
able to demonstrate its intrinsic behavior in Eq. (\ref{first:1})
through several oscillations.  Also, notice that the main peak of Eq.
(\ref{now}) vanishes at $\pi /\tau_{\rm max}$, so we should pick
\begin{equation}  
\frac{\pi}{T_{f}} \ll 
|\omega_0|  <
\frac{\pi}{\tau_{\rm max}} \label{C} 
\end{equation}
such that $\omega_0T_{f}\gg 1$, yet $\omega_0 \tau_{\rm max} < \pi$,
and so
\begin{equation} A \approx \frac{1}{T_{f}} \langle |f(\omega_0)|^2 \rangle
\end{equation}
by Eq. (\ref{puff}), (\ref{now}), where $\langle \rangle$ is taking
the proper average over $\omega_0$'s that satisfy Eq. (\ref{C}).  As
both $T_{f}$ and $\tau_{\rm max}$ goes to infinity, $\omega_0
\rightarrow 0$, which is the true meaning of ``taking $\omega=0$'' in
$|f(\omega)|^2/T_{f}$.

For equilibrium systems, $g(-\tau)= g^*(\tau)$. If the signal is
real, then $g(\tau)$ is also real, and we can write
\begin{equation} 
\int^{\tau_{\rm max}}_{0} g(\tau) d\tau 
\approx \frac{1}{2T_{f}} \langle |f(\omega_0)|^2 \rangle. 
\end{equation}

\section{Discrete Implementation}
\label{sec:discrete}

Fast Fourier Transform is the discrete implementation of Eq.
(\ref{first:1}) and (\ref{first:2}).  Let
\begin{equation} 
f_k = f\left( (k-1) \triangle\right), \; \triangle = \frac{T_{f}}{N}, \;  k=1..N,
\end{equation}
then 
\begin{eqnarray}
F_n =&& \sum^{N}_{k=1} f_k 
\exp\left(i\frac{2\pi (n-1)(k-1)}{N} \right ), \nonumber\\
f_k =&& \frac{1}{N} \sum^{N}_{n=1} F_n 
\exp\left(-i\frac{2\pi (n-1)(k-1)}{N}\right).
\end{eqnarray}
Obviously, $F_n$ corresponds to $\frac{1}{\triangle}f\left(\omega =
  \frac{2\pi (n-1)}{N\triangle}\right)$ when $\omega \triangle \ll 1$,
i.e., when the sample frequency is much bigger than the frequency of
interest. As as a convention, $F_0$ is always set to zero after the
Fast Fourier Transform by us, because it represents the heat current
average during the simulation, which should be subtracted off from the
signal.

Let us define
\begin{equation}
g_j = \frac{1}{N} \sum^{N}_{k=1} f^*_k f_{k+j-1}, \;\; j = 1..N, 
\end{equation}
where $f_{k+j-1 > N}$ is interpreted as $f_{k+j-1-N}$, and
\begin{eqnarray}
G_n &&= \sum^{N}_{j=1} g_j \exp\left( i\frac{2\pi (n-1)(j-1) }{N}
    \right ) \nonumber\\
    &&= \frac{1}{N} \sum^{N}_{j,k=1} f^*_k f_{k+j-1} 
\exp\left( i\frac{2\pi (n-1)(j-1) }{N}
    \right )
\nonumber\\
    &&= \frac{1}{N} |F_n|^2.
\end{eqnarray}
And so,
\begin{equation} 
g_j = \frac{1}{N} \sum^{N}_{n=1} G_n 
\exp\left(-i\frac{2\pi (n-1)(j-1)}{N}\right).  
\label{discrete_inverse}
\end{equation}
The idea is that we can calculate the correlation function by first
FFT the original signal $f_k$, take the square module of it and divide
by $N$, than inverse FFT back.

The discrete power spectrum summation
\begin{equation}
\frac{|F_n|^2}{N} = \frac{1}{N} \sum^{N}_{k,k^\prime=1} f_k f_{k^\prime}^*
\exp\left(i\frac{2\pi (n-1)(k^\prime-k)}{N}\right)
\label{corr}
\end{equation}
is represented graphically in Fig. \ref{noise}. The discretized
maximum correlation length is $K = N\tau_{\rm max}/T_f$. The shaded
region lies between $|k-k^\prime| < K$, so $f_k, f_{k^\prime}$ are
correlated within. Outside the shaded region their products can be
regarded as random noise, so passing through a wave with enough many
oscillations in the range would filter them out, which then requires
\begin{equation}
\frac{2\pi n}{N}N  = 2\pi n \gg 2\pi. \label{A}\end{equation}

We may also want to sum up the correlation $f_k f_{k^\prime}^*$ inside
the shaded region, which requires the phase factor $\exp\left(i2\pi
  (n-1)(k^\prime-k)/N\right)$ in Eq.  (\ref{corr}) to be approximately
unity for $|k-k^\prime|<K$. Thus if
\begin{equation}
\frac{2\pi n}{N}K \ll \frac{\pi}{2},
\label{B}
\end{equation}
or combined with Eq. (\ref{A}),
\begin{equation} 1\ll n \ll N/4K, \label{BB}
\end{equation}
we would have
\begin{equation}
\frac{|F_n|^2}{N} \approx 2 \langle \sum^{K}_{j=1} f_1 f_{j} \rangle.
\label{dudu}
\end{equation}
It is easy to see that finding such an $n$ is equivalent to choosing
$\omega_{0}$ in the case of continuous variables.

\section{Simulation results and correlation integral termination}
\label{sec:termination}

From now we will properly redefine ${\bf J}^q(t)$ such that the $1 /
k_{\rm B}T^2\Omega$ factor in Eq. (\ref{first}) vanishes.  Also,
$\tau$ will take the unit of $ps$ and $\underline{\underline{\kappa}}$
will take the unit of W/m/K. The spectral analysis of $f(t)$ can be
easily generalized to vector ${\bf J}^q(t)$ whereas we define
\begin{equation}
g_{ij}(\tau) = \langle J^q_i(0) J^q_j(\tau) \rangle,
\end{equation}
so there is
\begin{equation}
\kappa_{ij} = \int_{0}^{\infty}d\tau g_{ij}(\tau).
\end{equation}
${\bf J}^q(t)$ and $g_{ij}(\tau)$ are real, and $g_{ij}(\tau) =
g_{ji}(-\tau)$.  Furthermore, by time-reversal invariance of physical
systems, $g_{ij}(\tau) = g_{ij}(-\tau)$. So we should have
\begin{eqnarray} 
&& g_{ij}(\tau) \nonumber \\
=&& \frac{1}{2\pi T_{f}} \int^{+\infty}_{-\infty} 
d\omega \langle J_i^*(\omega) J_j(\omega) \rangle 
\exp (-i\omega \tau) \nonumber \\
=&& \frac{1}{2\pi T_{f}} \int^{+\infty}_{-\infty} 
d\omega {\rm Re} [\langle  J_i^*(\omega) J_j(\omega) \rangle]
\exp (-i\omega \tau).
\label{Gij}
\end{eqnarray}

For crystals with point symmetry higher than ${\rm T}_{\rm d}$, such
as $\beta$-SiC, $\underline{\underline{\kappa}}$ is a scalar matrix:
$\kappa_{ij} = \kappa \delta_{ij}$, and we can use this symmetry by
defining $|J(\omega)|^2 = (|J^q_{xx}(\omega)|^2 + |J^q_{yy}(\omega)|^2
+ |J^q_{zz}(\omega)|^2)/3$ and $g(\tau) = (g_{xx} + g_{yy} +
g_{zz})/3$, from the data of a single MD run.  So $\kappa =
(\kappa_{xx} + \kappa_{yy} + \kappa_{zz})/3$.

Theoretically the power spectrum $|J(\omega)|^2/2T_f$ provides complete
information about the thermal conductivity, because
\begin{eqnarray} 
\kappa = && \int_0^{+\infty} g(\tau) d\tau \nonumber\\
=&& \lim_{\omega_0\rightarrow 0} \frac{1}{2T_{f}} 
\left | J(\omega_0) \right |^2.
\label{AA}
\end{eqnarray}
In practice however, it is not a straightforward implementation
because $|J(\omega)|^2/2T_{f}$ from a single MD run is usually a
``cloud'' of points which requires smoothing. Fig.  \ref{freq5} shows
such a power spectrum from a 2 million timestep MD simulation of a
$\beta$-SiC crystal in a 512-particle cell at 1500K. The $x$-axis is
the frequency number $n$ after FFT. The $y$-axis is properly scaled to
the thermal conductivity unit W/m/K. The raw values are plotted in
empty circles, while results after properly smoothing are shown in
stars with solid line to guide the eye. As we show in the last
section, taking $\omega_0$ to zero in Eq.  (\ref{AA}) is equivalent to
choosing an $n$ on the discrete grid such that $1\ll n \ll N/4K$. We
shall see in Fig.  \ref{corr5} that the correlation $\langle
J(0)J(\tau) \rangle$ is significant for $\tau < \tau_{\rm max} \approx
20$ $ps$, which corresponds to $K = 60,000$.  So $n$ should satisfy $1
\ll n \ll 10$, say $n=6$ (in this case because $|J(\omega)|^2$ is the
average of three independent directions, the bounds are not so
strict), which roughly corresponds on Fig.  \ref{freq5} to a thermal
conductivity value of $70\pm 15$ W/m/K, in reasonably good agreement
with the experimental value of $62.6$ W/m/K.  However, this is
probably the most we can do from the power spectrum; the indeterminate
nature of this inference procedure forces us to find better, or at
least more automatic, ways of finding out the simulation thermal
conductivity result.

Notice that the requirement Eq. ($\ref{BB}$) is only barely satisfied
in the above example due to the smallness of $N/4K$, which is the
reason we call this a {\em minimal run}. As we shall demonstrate
later, the error bar of one's result is largely determined by the
ratio $N/4K$. By doing longer simulation runs at the same temperature,
one increases $N/4K$ and thus improves the quality of his result.

As a brief digression, we note (see inset graph of Fig. \ref{freq5})
that the power spectra sometimes show small peaks at finite
frequencies.  Besides roughening up $g(\tau)$ a bit, this will have
little influence on the calculated thermal conductivity.
Nevertheless, the nature of these peaks, whether to be something
physical or as a simulation artifact such as finite cell size effect,
remains unknown.

We now turn to the correlation function picture by inverse Fourier
transforming $|{J}(\omega)|^2$ using Eq.  (\ref{discrete_inverse}).
Fig. \ref{corr5} shows in solid line the normalized correlation
function $g(\tau)/g(0)$ after proper local smoothing.  Only the first
$1/8$ portion of entire data is plotted in the graph (total run length
= $747$ $ps$) because it appears to be completely random fluctuations
after that.  Clearly, we need a termination rule if somehow we want to
retrieve the $\kappa = \int_0^{+\infty} g(\tau) d\tau$ information.
Fig. \ref{corr5} represents the typical outcome of a minimal run,
which provides a test example for the termination criterion.

One observation is that in the case of $\beta$-SiC, especially at high
temperatures, the {\em actual} shape of $g(\tau)$ closely resembles
that of an exponential decay with somewhat faster relaxation in the
first few pico-seconds and a slight flattening of the tail. This
observation are confirmed by simulation results after multiple runs
are done and better statistics are collected (see Fig.
\ref{per_corr_1700}, \ref{per_corr_760}, \ref{per_corr_284}).  It
induces us to propose the following exponential fitting or EF
criterion: we can do a least square fit of $g(\tau)$ to an exponential
decay form $ae^{-(\tau-\tau_f)/b}$ in the time range $[\tau_i,
\tau_f]$, where simulation accuracy is still believed to be good, in
order to determine the best coefficients $\{a, b\}$.  We then
numerically integrate $g(\tau)$ up to $\tau_f$, and add with the
analytically estimated tail contribution $ab$. This termination rule
only requires two parameters to be chosen beforehand, $\tau_i$ and
$\tau_f$, depending on simulation run length, correlation time and
number of runs being averaged. The recommended values for $\beta$-SiC
are given in Table \ref{param}, along with other simulation
parameters.

Another more straightforward way is to postulate that the actual heat
current correlation function $g(\tau)$ {\em never} turns negative. We
know that certain self-correlation functions such as the velocity
autocorrelation function in liquids can have negative parts sometimes,
but the occurrence is rare and it is very unlikely that it happens in
solid SiC. If this is the case, then at the point where the calculated
correlation curve first turns negative, one can argue that the
calculated curve must be finally corrupted by random fluctuations, and
data from now on can not be trusted. One can then estimate the thermal
conductivity by numerically integrating the raw data up to the point
where the ``first dip'' happens, and this is called the first-dip or
FD rule.

Usually $g(\tau)$ from a long simulation run with a large $N/4K$ ratio
or the average of multiple runs will give FD and EF results that are
close to each other, especially at high temperatures; and large
difference occurs for those runs that are not well-converged (small
$N/4K$) and can be seen to be rather irregular. We list several
examples down below:

\begin{itemize}
  
\item 2 million step run, $T=1500$K, 512-particle cell (see Fig.
  \ref{freq5}, \ref{corr5}). EF result: 71.2 W/m/K; FD result: 67.8
  W/m/K; Exp't: 62.6 W/m/K.
  
\item 2 million step run, $T=1055$K, 216-particle cell.  EF
  result: 88.6 W/m/K; FD result: 89.3 W/m/K; Exp't: 95.7 W/m/K.

\item 2 million step run, $T=790$K, 216-particle cell. EF result:
  172.6 W/m/K; FD result: 229.7 W/m/K; Exp't:  162.4 W/m/K.
  
\item Average of $4\times 4$ million step runs (see Fig.
  \ref{per_corr_1700}),  $T=1700$K, 216-particle cell. EF result:
  69.2 W/m/K; FD result: 77.5 W/m/K; Exp't: 51.9 W/m/K.
  
\item Average of $4\times 4$ million step runs (see Fig.
  \ref{per_corr_760}),  $T=760$K, 216-particle cell. EF result:
  152.2 W/m/K; FD result: 162.7 W/m/K; Exp't: 138.5 W/m/K.
  
\item Average of $12\times 8$ million step runs (see Fig.
  \ref{per_corr_284}),  $T=284$K, 216-particle cell. EF result:
  146.2 W/m/K; FD result: 161.3 W/m/K; Exp't: 318.6 W/m/K.

\end{itemize}

Just by looking at the above reports and comparing with experiments,
one cannot clearly distinguish between the error due to the {\em
  model}, such as inexact interatomic potential or the negligence of
quantum effects which is quite damaging at low temperatures\cite{Li},
and the error due to the {\em implementation}, such as the termination
criteria. They can only be distinguished after doing many identical
runs with different random seeds, and accumulating enough data so that
we can claim certain result to be the {\em true} simulation result,
and then benchmark the quality of a specific implementation giving it
only part of the information. This procedure would be very 
time-consuming.  By experience, we think that for a single minimal
run, the EF result may be more accurate because it acts like a filter;
but when one can manage to do long simulation runs, and many of them,
the FD rule may yield better results because one then know the tail
part to be real. Table \ref{param} lists the parameters which we think
are appropriate for $\beta$-SiC crystal at various temperatures.

Finally, we note that the difficulty in benchmarking implementations
does not exist if we know the exact correlation function. This is done
in Section \ref{sec:model} where we propose a simple analytically
solvable model.


\section{Simple Solvable Model}
\label{sec:model}

Consider discrete series
\begin{equation}
a_0 = 0, \;\;\;\; a_i = a_{i-1}e^{-1/L}+ R_{i},
\label{sequence}
\end{equation}
where $R_i$ is an independent random number uniformly distributed on
$(-0.5, 0.5)$, and $L$ is an arbitrary constant signifying the
characteristic correlation length. The series can be shown to have a
simple correlation function:
\begin{eqnarray} 
g_k =&& \langle a_{i} a_{i+k-1} \rangle \nonumber\\
    =&& \langle a_{i} a_{i+k-2}e^{-1/L}+a_i R_{i+k-1}
    \rangle \nonumber\\
    =&& \langle a_{i}a_{i+k-2} \rangle e^{-1/L},
\end{eqnarray}
because $R_{i+k-1}$ is not correlated with $a_i$. And so
\begin{equation} 
g_k = \langle a_{i}a_{i} \rangle e^{-(k-1)/L}.
\end{equation}
Now,
\begin{eqnarray} 
a_{i} =&& R_i+a_{i-1}e^{-1/L} \nonumber\\
      =&& R_i+R_{i-1}e^{-1/L}+R_{i-2}e^{-2/L}+.. 
\end{eqnarray}
so
\begin{eqnarray} 
\langle a_{i} a_{i} \rangle 
=&& \langle R_i^2 \rangle+\langle R_{i-1}^2 \rangle e^{-2/L}
+\langle R_{i-2}^2 \rangle e^{-4/L}+..  \nonumber\\
=&& \frac{\langle R^2 \rangle}{1-e^{-\frac{2}{L}}}. 
\end{eqnarray} 
Because 
$$
\langle R^2 \rangle = \int_{-0.5}^{0.5}x^2 dx = \frac{1}{12}, $$
there is
\begin{equation} 
\langle a_{i}^2 \rangle = \frac{1}{12(1-e^{-\frac{2}{L}})},
\end{equation}
and
\begin{equation} 
g_k = \frac{e^{-(k-1)/L}}{12(1-e^{-\frac{2}{L}})}.
\label{model_corr}
\end{equation}

The theoretical thermal conductivity is then
\begin{equation} 
\kappa \equiv \sum_{k=1}^{+\infty} g_k = 
\frac{1}{12(1-e^{-\frac{2}{L}})(1-e^{-\frac{1}{L}})}.
\end{equation}

Using this known result we can check our codes by directly plugging in
the series, and see if after the procedure in Fig.  \ref{flowchart},
we indeed get a correlation function of the correct form and a correct
``thermal conductivity'' value. This is the best way to ensure that
all the coefficients are made right.

Let us illustrate the issues discussed in this paper by numerically
studying three examples of Eq. (\ref{sequence}), following the steps
outlined in Fig. \ref{flowchart}:
\begin{itemize}
\item Case A: $N = 2^{21}$, $N/L = 1000$, representing a well-converged
  MD run. The theoretical thermal conductivity is $\kappa_{A} = 1.834
  \times 10^5$.
\item Case B: $N = 2^{21}$, $N/L = 100$, representing a minimal MD
  run\cite{explain_minimal}.  The theoretical thermal conductivity
  should be $\kappa_{B} = 1.833 \times 10^7$.
\item Case C: Identical to case B but with a different random number
  seed.
\end{itemize}
where we take care to let the sequence iterate $N/8$ steps to
``equilibrate'' before collecting the $N$ data points.  Fig.
\ref{signal} shows a small portion of the signal in case A for $5L
\approx 10,000$ steps.

Fig. \ref{spectra1000} plots the power spectrum $|F_n|^2/2N$ of case A
(the inset shows the raw data on a larger scale) versus the frequency
number $n$ between $[1, N/L]$, which are the only relevant ones to
thermal conductivity as we demonstrated in Section
\ref{sec:discrete}\cite{explain_minimal}.  The empty circles are the
raw values, while the stars are the smoothed result after passing
through a zero-phase filter of nearby 10 points (Matlab {\em
  filtfilt}), with solid line to guide the eye. We can see that the
raw data is still a ``cloud'' of points, very much like that of Fig.
\ref{freq5}, but after smoothing it becomes more regular, with
discernible leveling off as $n \rightarrow 0$, which is related to the
fact that as $L\gg 1$ $|F_n|^2/2N$ becomes a Lorenzian. The $n
\rightarrow 0$ limit, as we can see, corresponds to a thermal
conductivity value of about $(1.7\pm 0.3)\times 10^5$ at frequency
number $n$ around $N/25L$. This verifies Eq. (\ref{dudu}) and Eq.
(\ref{BB}) and the agreement is actually pretty good, for this
``well-converged run'' with $N/L = 1000$.

We then inverse FFT the power spectrum to get $g_k$. The numerical
result is plotted in Fig. \ref{corr1000} in solid line, against the
theoretical prediction of Eq. (\ref{model_corr}) in dashed line.
Since the only timescale in this problem is $L$, the $x$-axis is
scaled as $k/L$. The numerical result is also shown in the inset on a
larger scale. We can clearly discern a region ($5<k/L<8$) in the
numerical result where $g_k \approx 0$ and is almost flat and parallel
to the $x$-axis, to the left of which $g_k$ is regular and
well-converged, and to the right it begins to randomly fluctuate.  Due
to the existence of this region, the thermal conductivity result we
get using $g_k$ (versus the power spectrum method) is bound to be very
good, no matter which termination rule is invoked. And this is what
actually happens in defect-state SiC thermal conductivity
calculations.

We then apply our EF rule to the data: the region of fitting
($[\tau_i, \tau_f]$ in last section) is chosen to be $k/L = [0.5,
1.5]$ (indicated by the two stars) which roughly corresponds to same
position in $\beta$-SiC. The resultant exponential curve is plotted in
dotted lines, which because it {\em has} the correct functional
dependence, is almost indistinguishable from the theoretical curve.
The inferred characteristic decay time is $L^\prime / L = 1.0043$, and
the inferred thermal conductivity is $\kappa_{A}^{\rm EF} / \kappa_A =
0.979$.

If we apply the FD rule to the same data, we find that the first dip
occurs at $k_{\rm dip} / L = 4.918$, which of course is completely
random, and the FD estimation gives a thermal conductivity value of
$\kappa_{A}^{\rm FD} / \kappa_A = 0.962$.

The entire procedure for case A registers $2.75 \times 10^8$ flops on
Matlab, while the brute force method requires computations on the
order of $NL = 4.40 \times 10^9$ flops, thus demonstrating the
efficiency of this approach. The comparison of cost will go even more
favorably for the spectral method in case B and C, which represents
more realistic scenarios.  Furthermore, the spectral method loses no
information about the original data; and one can gain extra knowledge
about the power spectrum.

Now let us do case B.  The meaning of symbols in the figures will be
the same as those in case A.  Fig. \ref{spectra100} shows the power
spectrum, where the stars are the smoothed result of nearby 5 points.
Because $N/L$ is small, we have much less relevant frequency numbers
than those in Fig. \ref{spectra100}. The $n\rightarrow 0$ limit gives
a thermal conductivity value of $(2 \pm 0.6) \times 10^7$ around
$n=7$, which agrees well with the theoretical value of $\kappa_{B} =
1.833 \times 10^7$, considering this is a minimal run.

Fig. \ref{corr100} shows the correlation function $g_k$ versus $k/L$.
We still fit the numerical result to an exponential form in the region
$0.5 \le k/L \le 1.5$. The inferred characteristic decay time is
$L^\prime / L = 1.280$, and the inferred thermal conductivity is
$\kappa_{B}^{\rm EF} / \kappa_B = 1.518$. The fitted exponential curve
is shown in dotted lines, where the region of fitting is between the
two stars. We can see that neither the original data nor the EF curve
agrees well with the theoretical curve, due to small $N/L$, i.e.,
insufficient data. However, the EF curve is able to filter out the
fictitious bulge between $4< k/L <8$.

If we apply the FD rule to the same data, we find that in this case
the first dip happens at $k_{\rm dip} / L = 8.385$, and the FD
estimation of the thermal conductivity value is $\kappa_{A}^{\rm FD} /
\kappa_A = 1.854$.

Neither of these results are good, but it is the kind of accuracy we
should expect from a single minimal run. The power spectrum estimation
in fact does better than the $g_k$ estimations in this case, but as we
said, the procedure depends more heavily on human, which is usually
not a good idea.

Now let us turn to case C, whose conditions are completely identical
with case B except changing the initial random number generator seed
to have a different sequence.  The goal is to find out how minimal
runs differ from case to case.  As it turns out, their difference can
be quite large. Fig. \ref{spectra100_c} shows the power spectrum of
case C, which does not seem very different from Fig.
\ref{spectra100}.  The $n\rightarrow 0$ limit gives a thermal
conductivity value of $(1.7 \pm 0.4) \times 10^7$ at about $n=8$.
However, the correlation function (Fig. \ref{corr100_c}) behaves much
better at small $k/L$ than case B. Applying the EF rule between $0.5
\le k/L \le 1.5$, we find that $L^\prime/L = 0.778$, and the inferred
thermal conductivity value is $\kappa_{C}^{\rm EF} / \kappa_C =
0.972$. If we apply the FD rule, the first-dip in correlation function
happens at $k_{\rm dip} / L = 5.786$, and the FD estimation gives
$\kappa_{A}^{\rm FD} / \kappa_A = 1.123$. So in this case both results
are better than the power spectrum estimation.  This might be due to
the fact the power spectrum in case C happens to fluctuate less at
small frequency numbers (compare inset graphs of Fig.
\ref{spectra100_c} with Fig.  \ref{spectra100}).

Now we had a glimpse of the subtleties involved in calculating the
thermal conductivity from minimal runs, which currently are often the
only things one can afford to do for realistic crystalline solids
($N/K \sim 30$ or $N/L \sim 100$). On the one hand, one expects the
error to be big, say $50\%$, and often biased to the positive side.
On the other hand, some minimal runs may happen to hit right on
target.  In general the difference from run to run is large.  The
difficulty is that unlike in this case where one knows the exact
answer, one cannot tell for sure whose error is smaller.  It is true
that certain correlation curves ``look'' better than others, such as
Fig.  \ref{corr100_c} than Fig.  \ref{corr100}, so one may assign
bigger weights to these results, but it is a complicated issue.

In the case of $\beta$-SiC, by making the assumption that the
correlation function follows an exponential decay form, or that it
never turns negative, we devise the EF and the FD termination rules.
As we said in Section \ref{sec:termination}, those curves which give
similar EF and FD values are often found to be closer to the correct
(well-averaged) result. This is also found to be the case for our
simple solvable model (compare case B and C); and it may also be
helpful for other system calculations.

Notice that the EF results in case A, B and C are always better than
the FD results, which we think is because the true correlation
function {\em is} exponentially decaying. In general, we believe that
the EF rule is more reliable for the data from a single minimal run,
while the FD rule is better for well-converged or well-averaged data.

Should people do minimal-run simulations? Sometimes they are forced
to: a minimal run for $\beta$-SiC at $T=500$K is 216 particles running
8 milllion steps under the Tersoff potential, which takes 50 processor
hours on the high-end computer of SGI Origin 2000.  The memory
requirement to analyze the data afterwards is also demanding.  But
even if one can do very long simulations, he should investigate before
doing so. As we saw, case A has a much better convergence than case B
and C, but it takes (relatively) 10 times longer to run.  We think
that if one does 10 minimal runs ($N=2^{21}/10$) and average over the
correlation functions, he will get comparable accuracy as case A. In
reality one may benefit more from the latter route because he has
taken a better ensemble average, by starting from different initial
conditions (see discussion in Section \ref{sec:introduction}).

\acknowledgments The development of the molecular dynamics code using
the Tersoff potential was supported in part by AFOSR Grant No.
F49620-96-1-0447.  Acknowledgment is also made to the donors of The
Petroleum Research Fund, administered by the ACS, for early support of
J. L., and to Honda Motor Company for the donation of the SGI Origin
2000 computer, on which parts of the calculation were carried out.

\begin{references}

\bibitem{Li} J. Li, L. J. Porter, S. Yip, J. Nucl. Mater. 246 (1998)
  53.

\bibitem{Green-Kubo} M. S. Green, J. Chem. Phys. 22 (1954) 398; R.
  Kubo, Rep. Prog. Phys. 29 (1986) 255.
  
\bibitem{Ladd} A. J. C. Ladd, W. Moran, W. G. Hoover, Phys. Rev. B 34
  (1986) 5058.
 
\bibitem{Lee} Y. H. Lee, R. Biswas, C. M. Soukoulis, C. Z. Wang, C. T.
  Chan, K. M. Ho, Phys. Rev. B 43 (1991) 6573.
 
\bibitem{Kitagawa} H. Kitagawa, Y. Shibutani, S. Ogata, Modell. Simul.
  Mater. Sci. Eng. 3 (1995) 521.
 
\bibitem{Kaburaki} H. Kaburaki, J. Li, S. Yip, to be published.
  
\bibitem{explain_g} \underline{\underline{g}}(0) can be thought of as
  a static thermodynamic susceptibility which is usually not as
  sensitively dependent on the material state.
 
\bibitem{Tersoff} J.  Tersoff, Phys. Rev. Lett. 64 (1990) 1757; J.
  Tersoff, Phys. Rev. B 49 (1994) 16349.
 
\bibitem{Press} W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P.
  Flannery,{\em Numerical Recipes in C -- the Art of Scientific
    Computing} (Cambridge University Press, Cambridge, 1992).
 
\bibitem{explain_minimal} There is some numerical difference between
  $K$ in Eq. (\ref{BB}) and $L$ in Eq. (\ref{sequence}): K signifies
  the {\em maximum} correlation length, above which correlation can be
  neglected, while $L$ is the characteristic period by which
  correlation decays by $e^{-1}$.  If we consider $e^{-3} \approx
  0.05$ to be ``negligibly small'', then $K \approx 3L$. For instance,
  in the $T=1500$K $\beta$-SiC example in Section
  \ref{sec:termination}, the fitted EF decay time is 5.95 $ps$ ($L$),
  while we consider correlation to be significant up to 20 $ps$ ($K$)
  by looking at Fig.  \ref{corr5}.  So there should be $1 \ll n \ll
  N/4K \approx N/12L$.
  
\end{references}

\widetext

\newpage
\begin{table}
\caption{Recommended Parameters for Perfect Crystal
  $\beta$-SiC Thermal Conductivity Calculation}
\begin{tabular}{|c|c|c|c|c|c|}
\hline
T (K) & Steps     & Timestep ($fs$) & Runs
      & EF $\tau_i$($ps$) & EF $\tau_f$($ps$) \\ \hline
250 -- 500           & $8\times 10^6$ &  0.30-0.35         &  6-10
      &  5.0             &  35.0          \\
500 -- 1000          & $4\times 10^6$ &  0.30-0.35         &  4-6
      &  3.5             &  21.0           \\
1000 -- 1500         & $4\times 10^6$ &  0.25-0.30         &  3-4
      &  1.0             &  9.0            \\
 $>$ 1500              & $2\times 10^6$ &   $<$ 0.25           &  3-4
      &  1.0             &  9.0            \\
\hline
\end{tabular}
\label{param}
\end{table}
\vspace{10mm} 

\newpage

\begin{figure}
\psB{flowchart.ps}
% \vspace{50mm} 
\caption{Flowchart of the spectral method to 
calculate the thermal conductivity.}
% \vspace{10mm}
\label{flowchart}
\end{figure}

\begin{figure}
\psB{noise.ps}
\caption{Discrete power spectrum summation.}
\label{noise}
\end{figure}

\begin{figure}
\psB{freq5.ps}
\caption{Power spectrum 
  (plotted in circles and properly normalized to the unit of W/m/K)
  versus frequency number $n$, of the heat current data from a 747
  $ps$ MD simulation of $\beta$-SiC crystal at $T=1500$K in a
  512-particle cell. The stars are results after 4-point smoothing,
  with solid line to guide the eye.}
\label{freq5}
\end{figure}

\begin{figure}
\psB{corr5.ps}
\caption{Normalized heat current correlation function $g(\tau)/g(0)$
  from a single MD run ($g(0) = 17.3$ W/m/K/$ps$), for perfect
  $\beta$-SiC at 1500K. The dashed line is the fit to exponential
  decay in the range 1 to 9 $ps$.}
\label{corr5}
\end{figure}

\begin{figure}
\psB{per_corr_1700.ps}
\caption{Well converged normalized heat current correlation function
  $g(\tau)/g(0)$ for perfect $\beta$-SiC at 1700K ($g(0) = 22.0$
  W/m/K/$ps$). The dashed line is the fit to exponential decay in the
  range 1 to 9 $ps$.}
\label{per_corr_1700}
\end{figure}

\begin{figure}
\psB{per_corr_760.ps}
\caption{Well converged normalized heat current correlation function
  $g(\tau)/g(0)$ for perfect $\beta$-SiC at 760K ($g(0) = 15.5$
  W/m/K/$ps$). The dashed line is the fit to exponential decay in the
  range 5 to 25 $ps$.}
\label{per_corr_760}
\end{figure}

\begin{figure}
\psB{per_corr_284.ps}
\caption{Well converged normalized heat current correlation function 
  $g(\tau)/g(0)$ for perfect $\beta$-SiC at 284K ($g(0) = 7.8$
  W/m/K/$ps$). The dashed line is the fit to exponential decay in the
  range 5 to 35 $ps$.}
\label{per_corr_284}
\end{figure}

\begin{figure}
\psB{signal.ps}
\caption{A small portion of the random signal $a_i$ in case A, with
  $N=2^{21}$ and $N/L=1000$.}
\label{signal}
\end{figure}

\begin{figure}
\psB{spectra1000.ps}
\caption{Power spectrum $|F_n|^2/2N$ versus frequency number $n$ in
  case A. The empty circles are the raw data, while the stars (with
  solid line to guide the eye) are the smoothed results after passing
  through a zero-phase filter of nearby 10 points. The inset graph
  shows the raw data on a larger scale.}
\label{spectra1000}
\end{figure}

\begin{figure}
\psB{corr1000.ps}
\caption{Correlation function $g_k$ versus $k/L$ in case A.
  The solid line is the numerical result; the dashed line is the
  theoretical prediction of Eq. (\protect\ref{model_corr}). The dotted
  line (indistinguishable from the theoretical curve in this case) is
  the EF fit to numerical result in the region $0.5 \le k/L \le 1.5$,
  indicated by the two stars. The inset graph shows the numerical
  result on a larger scale.}
\label{corr1000}
\end{figure}

\begin{figure}
\psB{spectra100.ps}
\caption{Power spectrum $|F_n|^2/2N$ versus frequency number $n$ in
  case B. The empty circles are the raw data, while the stars (with
  solid line to guide the eye) are the smoothed results after passing
  through a zero-phase filter of nearby 5 points. The inset graph
  shows the raw data on a larger scale.}
\label{spectra100}
\end{figure}

\begin{figure}
\psB{corr100.ps}
\caption{Correlation function $g_k$ versus $k/L$ in case B.
  The solid line is the numerical result; the dashed line is the
  theoretical prediction of Eq. (\protect\ref{model_corr}). The dotted
  line is the EF fit to numerical result in the region $0.5 \le k/L
  \le 1.5$, indicated by the two stars. The inset graph shows the
  numerical result on a larger scale.}
\label{corr100}
\end{figure}

\begin{figure}
\psB{spectra100_c.ps}
\caption{Power spectrum $|F_n|^2/2N$ versus frequency number $n$ in
  case C. The empty circles are the raw data, while the stars (with
  solid line to guide the eye) are the smoothed results after passing
  through a zero-phase filter of nearby 5 points. The inset graph
  shows the raw data on a larger scale.}
\label{spectra100_c}
\end{figure}

\begin{figure}
\psB{corr100_c.ps}
\caption{Correlation function $g_k$ versus $k/L$ in case C.
  The solid line is the numerical result; the dashed line is the
  theoretical prediction of Eq. (\protect\ref{model_corr}). The dotted
  line is the EF fit to numerical result in the region $0.5 \le k/L
  \le 1.5$, indicated by the two stars. The inset graph shows the
  numerical result on a larger scale.}
\label{corr100_c}
\end{figure}

\end{document}
